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Abstract 
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Q 

t^ \ Extended Bose Hubbard models with nearest neighbour interaction describe minimally the effect 



of long range interaction on ultra cold atoms in deep optical lattices. Rotation of such optical 

c/2 ■ 

C^ . lattices subject such neutral cold atoms to the effect of an artificial magnetic field. The modification 

of the phase boundaries of the density wave and Mott Insulator phases due to this rotation are 

C^ ' shown to be related to the edge spectrum of spinorial and scalar Harper equation. Corresponding 

^ ; 

"■ profiles of the checkerboard vortex states with sublattice modulated superfluid order parameter 



near density wave phase boundary are calculated. 
PACS numbers: 03.75.Lm, 64.70.Tg, 67.80.bd 



I. INTRODUCTION 

Ultra cold atomic condensates with short range interaction in deep optical lattices are 
described by the Bose Hubbard model [l| in the tight binding approximation and shows 
quantum phase transition from the superfiuid (SF) to Mott insulator (MI) phase due to the 
competition between nearest neighbour hopping and on site interaction 2|. Such condensates 
when subjected to an artificial magnetic field through rotation |3| or by imprinting motion 



fl 



dependent laser induced phases on their internal states [J], form vortices. The effect of an 
artificial magnetic field on the phases of cold atoms generated in the presence of an optical 
lattice |5| can either by studied by trapping more than one internal states of the atom in 
optical lattice |6|, 17| or by rotating the optical lattice [8[. This explores the effect of an 
artificial magnetic field on ultra cold neutral atoms in tight binding approximation. 

Extended Bose Hubbard (EBH) model 9|] that includes additionally interaction between 
atoms at different lattice sites described such cold atoms in optical lattices with long range 



interaction [10|]. Examples are dipolar cold atoms or polar molecules [ll|. In this paper we 
study the effect of rotation on such EBH model that includes nearest neighbour interaction 
(NNI) apart from the on-site interaction. The addition of the NNI to the Bose-Hubbard 
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hamiltonian has pronounced effect on the phases since the corresponding phase diagram 
Il4l | contains the density wave (DW) and supersolid (SS) phases apart from the MI and SF 
phases. The DW and MI phases lack coherence as the SF order parameter vanishes. Both 
have fixed number of particles at a given site. But DW has alternating particle numbers on 
successive sites ( Fig. [T] (a)) where as in the MI phase they are uniform. 

In the intriguing supersolid (SS) phase the superfiuid order parameter and the crystal 
order co-exist and the superfiuid density gets spatially modulated. The supersolid phase was 
first experimentally cited in solid helium [l5| though the interpretation of the experimental 
results was not without controversy [l6|. However if realized with cold atomic system in 
optical lattice, such a supersolid phase can be identified in a clearer fashion. An unambiguous 
way of identifying the SS phase is to study the modulation of the superfiuid order in the 
vortices created in such phase which will be different from the vortices created in an uniformly 
rotated superfiuid. To understand such vortex profiles one thus need to study the effect of 
such gauge field on the phases of EBH model. 

The phase diagram of ordinary BH model in presence of such gauge field or equivalently 



cold atoms in rotated optical lattice recently inspired a number of work 17H25 



he change 



in the nature of the quasiparticle excitations both near the phase boundary Il9l-|21L |24| ] as 
well as deep inside the superfluid phase due to the effect of the gauge field |25| has been 
studied extensively. However, the effect of gauge field on SS phases realized in EBH model 
received much less attention. In this paper we report the modification of the DW-SS phase 
boundary and the novel vortex profiles in SS phase near such phase boundary due to such 
gauge field. 

We calculate the modification of the DW phase boundary in the mean field approximation 
by using a reduced basis ansatz for the Gutzwiller variational wavefunction. The minimiza- 
tion of the energy functional very close to the DW phase boundary shows that the superfluid 



order parameter satisfies a spinorial Harper equation 



Consequently the phase bound 



ary can be determined from the edge of a Hofstadter butterfly (HB) spectrum 27|. In the 
resulting vortices, the spatial profile of the superfluid density shows a checkerboard like two 
sublattice modulation with a relative phase winding between the superfluid order parameter 
defined on each of these sublattices. We discuss their possible experimental detection. 



II. THEORETICAL FRAMEWORK 




FIG. 1: (a) alternating particle number in density wave phase (b) superfluid order parameter in 
super solid phase on the sites of A (red) and B(green) sublattices. 



We consider a square optical lattice in two spatial dimension rotated in the plane about 
z axis. The corresponding tight binding Hamiltonian in the co-rotating frame with onsite 



interaction and NNI is given by 

H = —ty ^{alaj exp{iipij) + h.c.) 

+ -'^ni{ni-l) + V^hinj - ii'^Ui (l) 

Here hopping amplitude t, NNI strength V and chemical potential fi are expressed in unit 
of the on site repulsion energy U. () implies that site index i,j on the two dimensional 
square lattice are the nearest neighbours and a|, cij, and Ui are bosonic creation, annihilation 
and number operators for the i-th site. We neglect the effect of an overall trap potential 
assuming that it is sufficiently shallow and gets neutralized by the centrifugal force that 
normally happens in the bulk of the system. 

The phase factor ipij = J ' dr.A{r) with the effective vector potential A{r) = {m/h)[fl x 
r) = -Kv^xy—yx) in the symmetric gauge. The resulting artificial magnetic field is 2Vtz where 
VL is the frequency of rotation. In Landau gauge A{r) = 2{m/K)VLxy which is more suitable 
for the experimental set-up in ref. ^. 

The quantity ^ = "^ = —^ § dr ■ V^ij gives the number of circulation quanta through 
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a unit cell in the square lattice and is gauge invariant. For the u = - [p and q are co-prime) 
as the boson hops around a unit cell in the square lattice it acquires a non trivial phase 
factor exp{—2'Kiv). To achieve a winding which is integer multiple of the 27r, the boson 
should therefore hop around q such unit cell leading to the formation of a magnetic unit cell 



281]. This in turn implies that if we denote the phase of the bosonic wave function by the 
direction of an arrow then as one goes around such magnetic unit cell, the arrow will rotate 
p times and the wavefunction will have —p vorticity in a magnetic unit cell. The same thing 
will happen even if we start from some other lattice than the square lattice as long as the 
number of flux quanta goes through the unit cell will remain v. Thus it imposes a topological 
constraint which does not depend on the local features such as the lattice potential.. 

The ground state of the Hamiltonian ([T]) can be found by variational minimization of 
{^\H\^) with a Gutzwiller wave function |\E') = Hi Xln. /nl'^i)- The variational parameters 
/^ are the amplitudes for the Fock state \ni) with n particles at site i. A detailed analysis of 
such variational mean fleld approaches that is generally used to study the many body states 



of the Bose-Hubbard hamiltonian is given in reference 29|]. For two dimensional lattice and 



< y < 4 and t = 0, the system will go through an alternating sequence of DW phase 



with no and no — 1 particles at successive sites, followed by a MI phase with no particles per 
site where no = 1, 2, 3, ■ ■ ■ . As t increases a SS phase appears before the DW state makes 
transition to a uniform superfluid phase. 

A. Phase boundary for the non rotating case 

The phase boundary of the DW phase can be determined analytically by obtaining the 
energy of the particle-hole type excitations using a reduced basis variational ansatz for the 
Gutzwiller wave function near the phase boundary. The DW phase consists of two sublattices 
A and B ( Fig. [1] (a)) having fixed no and no — 1 particles per site. Thus it is convenient 
to decompose |*) = (|*^))(|^^)). Here |^^) = Uf^^ili^"^) with \i/j'^) = En/n-'l^-iA) 
with /;^ = (5„,„„. Similarly |M/^) = 11,^/^1^*^) with |^'«) = E™/4"l^.s) with /^f = 
Sm^no-i- For the non rotating case of ^2 = 0, very close to the DW phase boundary only the 



neighbouring Fock states are populated 20|, l2]J]. Thus for all i and j 






We set (C-i,/^,C+i) = {eiA,V^-ei^-ei^,e,A) and (/^i, /™^ /^^i) 



{eiB, a/1 — ^iB ~ ^2_B' ^2_b) with variational parameters eiA,iB,£2A,2B all are ^ 1 to ensure 
the normalization condition of states \ip^^), {ip^'^). Also for the brevity of the notation we 
have written n^^ as n and rriig as m in the above expressions. Minimization of the en- 
ergy with respect to these four parameters gives four equations. Their non trivial solution 
demands 
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n- 1) - ft + AVm 




{n- fi + AVm) 


Atyjnm 


At^n{m + 1) 





-At^Jm{n + 1) 

(m - 1) - /i + AVn 




-4tv/(m + l)(n + l) 


At^Jn■m 


Atsjmin + 1) 
-At^{n + l){m + l) 





-At^Jn{m + 1) 


{m- 11 + AVn) 



(2) 
A particle (^) or hole {^) like excitation from either site of A and B are respectively given 

hj Sp = n + AVm,e^ = -[{n-l)+AVm\, e^ = m + AVn,e^ = -[{m-l)+AVn]. Defining 



£p'h = ^p,'h =F yu Eq. ([2]) gives the relation [l3[ 

epfe^e^ - m' [{n + l)et + <] [{m + 1).^ + me^] = (3) 

The above equation determines the chemical potential /i at each t for a given strength 
V of the nearest neighbour interaction and will give the phase boundary. To understand 
significance of this equation in a better way we compare it with the similar results obtained 
within the framework of other mean field approaches such as time dependent Gutzwiller 



mean field theory 13|, |29| where also minimal perturbation around a perfect density wave 
state is considered in the Fock space basis. Now, when such particle or hole like excitation is 
created over a perfect density wave state, they do not remain localized at a site, but moves 
around the lattice to create a Bloch wave to minimize their energy. The kinetic energy of 
such a Bloch wave is given by e{k) = 2t{coskx + cos ky), where kr^, ky are the components of 
the Bloch wave vector. The excitation spectrum of such particle-hole like excitations with 
inite wavevector can be obtained within the time dependent Gutzwiller mean field theory 



13| as 

e^e^ete^ - e{kf [{n + l)et + ne^] [{m + l)e^ + mef ] = 

where k = {kx,ky}. The density wave boundary can again be retrieved by taking zero wave 
vector limit, namely k^ — )■ 0, /cy — )■ 0, e{k) = At, which expectedly reproduces our result in 
Eq. Q. We again emphasize that all the above displayed relations are for two dimensional 
square lattice, but can be generalized in other dimensions. 

In the next subsection we shall extend the above treatment for the rotating case and 
will show that the limiting particle-hole excitation spectrum that determines such phase 
boundary in presence of the finite rotation ( or magnetic field) is actually the edge of a 
Hofstadter butterfly like energy spectrum. 



B. Rotated case 



For the rotated case, fi 7^ 0, we have 



^r<p'i] 



{iA,iB) 
i=N i=N 



2' 

1=1 rii i=l rii 



{iA,iB) "A rriB 

The first, second and fourth term respectively gives the mean kinetic, on site and nearest 
neighbour energy and the summation over i in the second and third term includes both 
the sublattices. In all further description again for brevity n^ and rriB will be written 
as n and m The superfiuid order parameter on two sublattices (Fig. [1] (b)) are given by 
0^,0^ = {di^), (dig) whereas the DW order parameter is given by (— l)*[(^i) — ;^(Xli'^i)] ^n 
any site i on either sublattices. Near the DW phase boundary again only the neighbouring 



Fock states will get occupied. The corresponding variational parameters {fn-iyfk^yfn^ 



+1^ 



for u sites are [A^^A^^^^*, Jl - \A(j)'^\^{\XY\^ + |A^-"|2), A^^A0^^], and, for ib sites we write 



ifLUJL'JL^i) as [6l-A<PT, a/1 - |A</.'#|2(|<51-|2 + \5l-\^),6l-A<P'i]. The superfiuid or- 



der parameters on the two sublattices are respectively given by 0^ = ^„ -\/rH-T/^'**/^+i 



and (f)^B = Xlm v^"^ + 1/nf */m+i- From these definitions it can be shown 0^ = A0^^ + 
0((A(/.*/)3) with Xy^ = ^(1 - v^Ai^) and similarly (f)'i = A(f)'^ + 0{{A(j)'^f) with 
^2^ — 7^+1(1 ~ v^'^i'^) Thus if we neglect third and higher order corrections, A0a,b can be 
replaced by the superfiuid order parameter (I)a,b on the two sublattices. Substituting these 
replacements and the expressions for variational parameters in the Eq. IHwe obtain 



($1^1*) = -2tRe Y, [e''''^'^<PT<P'i] + Y^ ^"^ ^tl^^""^ [^ ~ 2v^|Al^| - \X^^\^] + |Ai^p |0*/|2 

{iA,iB) *A 

[m — fi + 4Vn) 



n + 1 



+ E 



m 



[1- 27^151^1 -|5r|2] +|5[ 



s|2 



6-p 



Eg (5) 



with i^G is the energy of the pure density wave state. To obtain the ground state energy, 
{'^\H\'^) is first minimized with respect to A^"* and 5^^ yielding A^"* = -y/n ""^^^^^ ; 6\^ = 

m—fi+4Vn 
l+li-4Vn 



^ T+^-4v" - Substituting the above expressions in (\E'|if|\i/) , and, setting 0^ = ^Jei^^ 
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FIG. 2: (a) Hofstadter butterfly: the energy (e) spectrum for tlie Eq. ([8]) for various i^{0, 1}. 
The upper edge (marked red) gives the boundary of the density wave and the Mott Insulator lobe 
as explained in the text, (b) The first DW and MI lobe as a function of t, fi, v in mean field 
approximation where t and [i are in the unit of U . (c) Cross section of the plots in (b) that shows 
the modification of the first two density wave lobe and the first Mott lobe at various values of 
circulation quanta v. In all these plots V has been taken as 0.2 in the unit of U . 



and 0^ 



e20B ^^'^ ^ 



where ei 



(n— /^+4ym) 
n+1 



n 



n—n+AVm 



n,n ^ m), gives us the energy functional £ near the DW phase boundary as 

T 



and e2 = ei(m — ?> 



^T 



IB* 

B 



£ = -f E 

{iA,iB) 



[n ■ cr] 



^'t 4>'i 



^-P 



_ Eg (6) 

lA is 

The unit vector n = cosipi^igX + sincyjj^j^y and cr = axX + ayij, where ax,y are the Pauli 
matrices. The reduced basis ansatz assumes very low superfiuid density (0a,b <^ !)• This is 
valid very close to the phase boundary. Thus £ contain terms only linear in the superfiuid 
density. This is unlike the Gross-Pitaevskii energy functional which contains terms quadratic 
in the superfiuid density and is valid deep inside the superfiuid regime. 

Minimization of the above energy functional with respect to 0^*, 0^* gives equations for 
the superfiuid order parameter that can be written as a spinorial Harper equation, 

T \ V iT 



{iA,iB) 



n ■ <j] 



^a" <P'E 



t 



^'t <^'i 



(7) 



^M'l 



Its solution can be written as (t>{x,y) ® exp(— z 
the following symmetric gauge Harper equation 26 1 



exp[i^- 



wliere 0(x, y) satisfies 



0(x + 1, y)e''"'y + (j){x - 1, y)e-^^"^^ 
+ 4>ix, y + l)e^'"'^" + 4>{x, y - l)e''^"" 



■~(j){x,y) 



(8) 



i in the right hand side of the Eq. (|8]) can be mapped on the eigenvalues e of HB 27 1 
spectrum plotted in Figl2](a). 



III. RESULTS AND DISCUSSION 



The edge of the HB spectrum (marked red in Fig. [2] (a)) gives the highest eigenvalue of 
the Eq. ([H]) as function of u{0, 1}. This corresponds to the minimum value of t = tc = /"^ 
with non vanishing (p for each given value of /i, and, hence the boundary of the DW phase at 
that particular u ( marked red in Fig. [2] (b)). Same observation holds true for MI boundary 
at the MI-SF transition in a rotating optical lattice 20|, |2l| . Setting m = nin the preceeding 
analysis the MI-SF transition in rotated lattice can be studied for Extended Bose Hubbard 
model. The phase boundary of the ordinary Bose Hubbard model under rotation or magnetic 
field can be retrieved by setting V = and also putting n = m in t he p receding analysis. 
The results obtained in this way matches with those in reference 20|, |2l||. As compared to 



the modification of the phase boundary of a MI phase in ordinary BH model here also the 
phase boundary of the DW phase extends as the strength of the gauge field u is enhanced. 
This is due to the stronger localization of the bosonic states by the increasing strength of 
the gauge field. However the superfluid order parameter of the excitations at the boundary 
of the DW phase are different from those near the MI boundary as we shall see in Fig|31 
The modification of DW as well as MI phase boundaries are plotted in Fig. [2](b) and (c). 
However it is important to note that the Fig. [2](c) only provides analytically the phase 
boundary of DW as well as MI phases within reduced basis ansatz and does not provide the 
phases themselves over the entire t — fi plane for various z/ unlike in the references 12H14 1 . 



Particularly in mean field approximation it can be obtained numerically by using the full 
Gutzwiller wavefunction as was done in [13 1. 

At t = tc and u = jjj each magnetic unit cell that consists oi L x L lattice sites, contains 
one single vortex of unit winding. The strong sublattice modulation of the superfluid density 
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around the vortex core is shown in Fig. [3](a) for L = 16. The DW order parameter given 
in Fig. ^h) becomes 1 at the vortex core and co exists alongside the superfluid order in 
the bulk. Since at t = tc the systems become a supersolid in the mean-field Gutzwiller 



approximation 13|, the vortex structure in Fig. ^a) correspond to the vortex structure just 
at this transition boundary. 

We know that in a Hofstadter butterfiy problem, for z/ = -, a given degenerate Landau 
level is broken into q bands for p fiuxes through a given magnetic unit cell. In our present 
case we have taken ^ = 256- "^^^ highest of these energy levels correspond to the critical value 
of t = tc at the phase boundary. Thus one may think that the eigenfunction for the lower 
energy levels that correspond to higher values of t can be related with the superfiuid phases 
away from the phase boundary of the DW state inside the supersolid regime. However this 
simplistic argument is not completely correct since the entire derivation presented above is 
only within the reduced basis ansatz, which is valid for t ~ t^- Nevertheless we also plot the 
eigenfunction corresponding to a band which is very close to the highest band in in Fig. |3]^c). 
This approximately depicts the superfiuid order parameter in a rotated supersolid phase for 
t > tc , but still very close to the DW phase boundary. This state, corresponding to the 
lower band of the same spectrum contains multiple vortices in a given magnetic unit cell and 
the winding number of these vortices could also be integers > 1. Such a vortex structure is 
plotted in Fig. ^c). For calculating vortex structure at higher t that corresponds to deep 
inside the supersolid phase, one needs to go beyond the reduced basis ansatz and includes 
the non linear terms due to superfiuid interaction. 

Experimental detection of such vortices near the phase boundary is possible with the 
presently available techniques. The sublattice modulation of the superfiuid density can be 
detected through the time of fiight measurement and studying the resulting interference 
pattern j^. To measure the detailed vortex structure in a magnetic unit cell one can use 



Bragg scattering technique [30|] which is sensitive to the spatial phase distribution of the 
initial state (Sli, direction of rotation [32| and thus provides us a robust signature of the 
vortex state. 

To conclude we showed that the modification of phase boundaries of an EBH model 
due to rotation induced artificial magnetic field can be derived from the edge spectrum of 
a spinorial Harper equation. From the spectrum of the same equation we have explicitly 
demonstrated within mean field theory how the superfiuid and crystal order co-exist in the 
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FIG. 3: (a) checkerboard vortices at the density wave ( |2, 1, 2, 1, • • • )) phase boundary (t = tc) 
corrrsponding to the highest eigenvalue (the edge) of the hofstadter butterfly spectrum for i' = 
Yg— fg . The direction of the arrow gives ipi^ ^ig where as the color axis gives the superfluid density. 
The superfluid density is normalized by the maximum superfluid density at the boundary, (b) 
corresponding DW order parameter (c) More complicated vortex structure corresponding to the 
higher value of t corresponding to a lower eigenvalue ((254(16 x 16 — 2)th band)( d) corresponding 
DW order parameter. 

vortex profile of a supersolid around a DW vortex core. This can be used to identify the 
exotic supersohd phase in cold atom experiments. The above calculation can be generalized 
for the other variants of the EBH model such as one that includes next nearest neighbour 
interaction and can motivate further study for such vortices by going beyond the mean field 
approximation. 

We thank G. V. Pai, D. Goldbaum, E. Mueller, K. Seshadri, J. Avron, O. Gat and E. 
Altman, S. Sinha for helpful correspondences. The work of RS is supported by CSIR, India 
and the work of SG is supported by the Planning unit of IIT Delhi. 
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